# type2-stationary.R
# Philips, Andrew Q. "How to Avoid Incorrect Inference (While Gaining Correct Ones) in Dynamic Models". Forthcoming at Political Science Research and Methods. 
# andrew.philips@colorado.edu
# 2/9/21
#
# Figures produced in this R script:
#   --SI, Figure 21: "yi0-xi0-stationary-srboxplot.pdf"
#   --SI, Figure 22: "yi0-xi0-stationary-lrboxplot.pdf"
#   --Main document, Figure 6: "yi0-xi0-stationary-rejecttruth.pdf"
#   --SI, Figure 16: "yi0-xi0-stationary-rejecttruth-evar5.pdf"
#   --SI, Figure 17: "yi0-xi0-stationary-reject0.pdf"
#   --SI, Figure 18: "yi0-xi0-stationary-reject0-evar5.pdf"
#   --SI, Figure 19: "yi0-xi0-stationary-mse-sr.pdf"
#   --SI, Figure 20: "yi0-xi0-stationary-mse-lr.pdf"
# -----------------------------------

rm(list = ls())
setwd("~/Dropbox/My_Folder/Papers-Projects/Wlezien-Enns-PSRM/Final-Feb2021/scripts/")
# -----------------------------------
library(dynamac) # for lagging
library(nlWaldTest)
library(ggplot2)
library(grid)
library(gridExtra)
library(dplyr)


set.seed(39201848)
# -------------------------------------------------------------
# EXPERIMENT 5: Y ~ I(0), X ~ I(0), and related
sims = 2000 # no. sims
N <- c(50, 250, 1000) # number of obs
e.var <- c(1, 5) # error variance in Y
B.x.true = 2
B.lx.true = c(1,-1)
phi.true = c(0.2, 0.8)
# LRMs are:
(B.x.true + B.lx.true[1])/(1 - phi.true[1]) # 3.75
(B.x.true + B.lx.true[2])/(1 - phi.true[1]) # 1.25
(B.x.true + B.lx.true[1])/(1 - phi.true[2]) # 15
(B.x.true + B.lx.true[2])/(1 - phi.true[2]) # 5

container.ardl <- matrix(NA, nrow = sims*length(N)*length(e.var)*length(phi.true)*length(B.lx.true), ncol = 11) # null matrix to hold output
# cols = sims, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var, phi.true, B.lx.true
container.ecm <- matrix(NA, nrow = sims*length(N)*length(e.var)*length(phi.true)*length(B.lx.true), ncol = 11) # null matrix to hold output
# cols = sims, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var, phi.true, B.lx.true
container.ldv <- matrix(NA, nrow = sims*length(N)*length(e.var)*length(phi.true)*length(B.lx.true), ncol = 11) # null matrix to hold output
# cols = sims, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var, phi.true, B.lx.true
container.static <- matrix(NA, nrow = sims*length(N)*length(e.var)*length(phi.true)*length(B.lx.true), ncol = 8) # null matrix to hold output
# cols = sims, N, Beta.x, UL.beta.x, LL.beta.x, e.var, phi.true, B.lx.true

row <- 1
prog <- txtProgressBar(min = 0, max = sims, style = 3)
for (s in 1:sims) {
  for (n in N) {
    for (ev in e.var) {
      for (p in phi.true) {
        for (b in B.lx.true) {
          setTxtProgressBar(prog, value = s)
          
          x <- rnorm(n + 100, sd = 1)
          y.error <- rnorm(n + 100, sd = sqrt(ev))
          y <- y.error
          for(i in 2:length(x)) {
            y[i] <- p*y[i-1] + B.x.true*x[i] + b*x[i-1] + y.error[i]
          }
          
          # ARDL model:
          res.ardl <- lm(y[101:length(y)] ~ lshift(y, 1)[101:length(y)] + x[101:length(x)] + lshift(x, 1)[101:length(x)])
          # ECM model:
          res.ecm <- lm(dshift(y)[101:length(y)] ~ lshift(y, 1)[101:length(y)] + dshift(x)[101:length(x)] + lshift(x, 1)[101:length(x)])
          # LDV model
          res.ldv <- lm(y[101:length(y)] ~ lshift(y, 1)[101:length(y)] + x[101:length(x)])
          # static model
          res.static <- lm(y[101:length(y)] ~ x[101:length(x)])
          
          # grab QOI (FORALL): ----
          container.ardl[row, 1] <- container.ecm[row, 1] <- container.ldv[row, 1] <- container.static[row, 1] <- s # no sims
          container.ardl[row, 2] <- container.ecm[row, 2] <- container.ldv[row, 2] <- container.static[row, 2] <- n # T
          container.ardl[row, 9] <- container.ecm[row, 9] <- container.ldv[row, 9] <- container.static[row, 6] <- ev # error variance of Y
          container.ardl[row, 10] <- container.ecm[row, 10] <- container.ldv[row, 10] <- container.static[row, 7] <- p # true value of phi
          container.ardl[row, 11] <- container.ecm[row, 11] <- container.ldv[row, 11] <- container.static[row, 8] <- b # true value of effect of l.x
          
          # QOI (ARDL)
          container.ardl[row, 3] <- coef(res.ardl)[3] # beta.x
          container.ardl[row, 4] <- confint.default(res.ardl, parm = c(3), level = 0.95)[2] # UL
          container.ardl[row, 5] <- confint.default(res.ardl, parm = c(3), level = 0.95)[1] # LL
          # LR effects (ARDL)
          lrm <- nlConfint(res.ardl, "(a[3] + a[4])/(1-a[2])", level = 0.95)
          container.ardl[row, 6] <- lrm[1] # LRM
          container.ardl[row, 7] <- lrm[3] # LRM UL
          container.ardl[row, 8] <- lrm[2] # LRM LL
          
          # QOI (ECM)
          container.ecm[row, 3] <- coef(res.ecm)[3] # beta.d.x
          container.ecm[row, 4] <- confint.default(res.ecm, parm = c(3), level = 0.95)[2] # UL
          container.ecm[row, 5] <- confint.default(res.ecm, parm = c(3), level = 0.95)[1] # LL
          # LR effects (ECM)
          lrm <- nlConfint(res.ecm, "(a[4])/(-a[2])", level = 0.95)
          container.ecm[row, 6] <- lrm[1] # LRM
          container.ecm[row, 7] <- lrm[3] # LRM UL
          container.ecm[row, 8] <- lrm[2] # LRM LL
          
          # QOI (LDV)
          container.ldv[row, 3] <- coef(res.ldv)[3]
          container.ldv[row, 4] <- confint.default(res.ldv, parm = c(3), level = 0.95)[2] # UL
          container.ldv[row, 5] <- confint.default(res.ldv, parm = c(3), level = 0.95)[1] # LL
          # LR effects (LDV)
          lrm <- nlConfint(res.ldv, "(a[3])/(1-a[2])", level = 0.95)
          container.ldv[row, 6] <- lrm[1] # LRM
          container.ldv[row, 7] <- lrm[3] # LRM UL
          container.ldv[row, 8] <- lrm[2] # LRM LL
          
          # QOI (static)
          container.static[row, 3] <- coef(res.static)[2]
          container.static[row, 4] <- confint.default(res.static, parm = c(2), level = 0.95)[2] # UL
          container.static[row, 5] <- confint.default(res.static, parm = c(2), level = 0.95)[1] # LL
          
          row <- row + 1
        }
      }
    }
  }
}

# Create rejection rates and label:
container.ardl <- as.data.frame(container.ardl)
# cols = sim, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var, phi, B.lx.true
colnames(container.ardl) <- c("sims", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var", "phi", "B.lx.true")
container.ecm <- as.data.frame(container.ecm)
# cols = sim, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var, phi, B.lx.true
colnames(container.ecm) <- c("sims", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var", "phi", "B.lx.true")
container.ldv <- as.data.frame(container.ldv)
# cols = sim, N, Beta.x, UL.beta.x, LL.beta.x, LR.x, LR.x.UL, LR.x.LL, e.var, phi, B.lx.true
colnames(container.ldv) <- c("sims", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "LR.x", "LR.x.UL", "LR.x.LL", "e.var", "phi", "B.lx.true")
container.static <- as.data.frame(container.static)
# cols = sim, N, Beta.x, UL.beta.x, e.var, phi, B.lx.true
colnames(container.static) <- c("sims", "Time", "Beta.x", "Beta.x.ul", "Beta.x.ll", "e.var", "phi", "B.lx.true")

# proportion correct rejection, SR (from 0)
container.ardl$reject.beta.x <- ifelse(container.ardl$Beta.x.ll > 0 | container.ardl$Beta.x.ul < 0, 1, 0)
container.ecm$reject.beta.x <- ifelse(container.ecm$Beta.x.ll > 0 | container.ecm$Beta.x.ul < 0, 1, 0)
container.ldv$reject.beta.x <- ifelse(container.ldv$Beta.x.ll > 0 | container.ldv$Beta.x.ul < 0, 1, 0)
container.static$reject.beta.x <- ifelse(container.static$Beta.x.ll > 0 | container.static$Beta.x.ul < 0, 1, 0)

# proportion rejection, SR (from truth)
container.ardl$reject.beta.x.truth <- ifelse(container.ardl$Beta.x.ll > B.x.true | container.ardl$Beta.x.ul < B.x.true, 1, 0)
container.ecm$reject.beta.x.truth <- ifelse(container.ecm$Beta.x.ll > B.x.true | container.ecm$Beta.x.ul < B.x.true, 1, 0)
container.ldv$reject.beta.x.truth <- ifelse(container.ldv$Beta.x.ll > B.x.true | container.ldv$Beta.x.ul < B.x.true, 1, 0)
container.static$reject.beta.x.truth <- ifelse(container.static$Beta.x.ll > B.x.true | container.static$Beta.x.ul < B.x.true, 1, 0)

# proportion correct rejection of 0, LR:
container.ardl$reject.lrm <- ifelse(container.ardl$LR.x.LL > 0 | container.ardl$LR.x.UL < 0, 1, 0)
container.ecm$reject.lrm <- ifelse(container.ecm$LR.x.LL > 0 | container.ecm$LR.x.UL < 0, 1, 0)
container.ldv$reject.lrm <- ifelse(container.ldv$LR.x.LL > 0 | container.ldv$LR.x.UL < 0, 1, 0)

# proportion correct rejection, LR (from truth):
container.ardl$reject.lrm.truth <- ifelse(container.ardl$LR.x.LL > ((B.x.true + container.ardl$B.lx.true)/(1-container.ardl$phi)) | container.ardl$LR.x.UL < ((B.x.true + container.ardl$B.lx.true)/(1-container.ardl$phi)), 1, 0)
container.ecm$reject.lrm.truth <- ifelse(container.ecm$LR.x.LL > ((B.x.true + container.ecm$B.lx.true)/(1-container.ecm$phi)) | container.ecm$LR.x.UL < ((B.x.true + container.ecm$B.lx.true)/(1-container.ecm$phi)), 1, 0)
container.ldv$reject.lrm.truth <- ifelse(container.ldv$LR.x.LL > ((B.x.true + container.ldv$B.lx.true)/(1-container.ldv$phi)) | container.ldv$LR.x.UL < ((B.x.true + container.ldv$B.lx.true)/(1-container.ldv$phi)), 1, 0)

# squared error:
container.ardl$sqerr.beta <- (container.ardl$Beta.x - B.x.true)^2
container.ecm$sqerr.beta <- (container.ecm$Beta.x - B.x.true)^2
container.ldv$sqerr.beta <- (container.ldv$Beta.x - B.x.true)^2
container.static$sqerr.beta <- (container.static$Beta.x - B.x.true)^2
container.ardl$sqerr.lrm <- (container.ardl$LR.x - ((B.x.true + container.ardl$B.lx.true)/(1-container.ardl$phi)))^2
container.ecm$sqerr.lrm <- (container.ecm$LR.x - ((B.x.true + container.ecm$B.lx.true)/(1-container.ecm$phi)))^2
container.ldv$sqerr.lrm <- (container.ldv$LR.x - ((B.x.true + container.ldv$B.lx.true)/(1-container.ldv$phi)))^2
table(container.ardl$reject.lrm.truth)
table(container.ecm$reject.lrm.truth)
table(container.ldv$reject.lrm.truth)

# collapse everything down:
container.collapse.ardl <- container.ardl %>%
  group_by(Time, e.var, phi, B.lx.true) %>%
  summarize(reject.beta.x = mean(reject.beta.x),
            reject.lrm = mean(reject.lrm),
            reject.beta.x.truth = mean(reject.beta.x.truth),
            reject.lrm.truth = mean(reject.lrm.truth),
            mse.beta = mean(sqerr.beta),
            mse.lrm = median(sqerr.lrm))

container.collapse.ecm <- container.ecm %>%
  group_by(Time, e.var, phi, B.lx.true) %>%
  summarize(reject.beta.x = mean(reject.beta.x),
            reject.lrm = mean(reject.lrm),
            reject.beta.x.truth = mean(reject.beta.x.truth),
            reject.lrm.truth = mean(reject.lrm.truth),
            mse.beta = mean(sqerr.beta),
            mse.lrm = median(sqerr.lrm))

container.collapse.ldv <- container.ldv %>%
  group_by(Time, e.var, phi, B.lx.true) %>%
  summarize(reject.beta.x = mean(reject.beta.x),
            reject.lrm = mean(reject.lrm),
            reject.beta.x.truth = mean(reject.beta.x.truth),
            reject.lrm.truth = mean(reject.lrm.truth),
            mse.beta = mean(sqerr.beta),
            mse.lrm = median(sqerr.lrm))

container.collapse.static <- container.static %>%
  group_by(Time, e.var, phi, B.lx.true) %>%
  summarize(reject.beta.x = mean(reject.beta.x),
            reject.beta.x.truth = mean(reject.beta.x.truth),
            mse.beta = mean(sqerr.beta))

#save.image("scenario-yi0-xi0-stationary-data.RData")
#load("scenario-yi0-xi0-stationary-data.RData")

# box plots of non-collapsed data, for the SI:
container.ardl$phi.label[container.ardl$phi == .8] <- container.ecm$phi.label[container.ecm$phi == .8] <- container.ldv$phi.label[container.ldv$phi == .8] <- container.static$phi.label[container.static$phi == .8] <- "Alpha = 0.8"
container.ardl$phi.label[container.ardl$phi == .2] <- container.ecm$phi.label[container.ecm$phi == .2] <- container.ldv$phi.label[container.ldv$phi == .2] <- container.static$phi.label[container.static$phi == .2] <- "Alpha = 0.2"

p1 <- ggplot(subset(container.ardl, e.var == 1), aes(x = as.factor(B.lx.true), y = Beta.x, fill = as.factor(Time) )) + geom_boxplot() + facet_wrap(~phi.label, scale = 'free') + labs(fill = "Time") + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = B.x.true, color = "black", size = 1) + ylab("Estimate") + xlab("Coefficient on lag of X") + ggtitle("ARDL")
p2 <- ggplot(subset(container.ecm, e.var == 1), aes(x = as.factor(B.lx.true), y = Beta.x, fill = as.factor(Time) )) + geom_boxplot() + facet_wrap(~phi.label, scale = 'free') + labs(fill = "Time") + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = B.x.true, color = "black", size = 1) + ylab("Estimate") + xlab("Coefficient on lag of X") + ggtitle("ECM")
p3 <- ggplot(subset(container.ldv, e.var == 1), aes(x = as.factor(B.lx.true), y = Beta.x, fill = as.factor(Time) )) + geom_boxplot() + facet_wrap(~phi.label, scale = 'free') + labs(fill = "Time") + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = B.x.true, color = "black", size = 1) + ylab("Estimate") + xlab("Coefficient on lag of X") + ggtitle("LDV")
p4 <- ggplot(subset(container.static, e.var == 1), aes(x = as.factor(B.lx.true), y = Beta.x, fill = as.factor(Time) )) + geom_boxplot() + facet_wrap(~phi.label, scale = 'free') + labs(fill = "Time") + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = B.x.true, color = "black", size = 1) + ylab("Estimate") + xlab("Coefficient on lag of X") + ggtitle("Static")
pdf("yi0-xi0-stationary-srboxplot.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, p3, p4, ncol = 2)
dev.off()

# and for the long-run. For this, center around true lR effect since that differs based on alpha and l.x coef:
container.ardl$centered.lrm <- container.ardl$LR.x - ((B.x.true + container.ardl$B.lx.true)/(1-container.ardl$phi))
container.ecm$centered.lrm <- container.ecm$LR.x - ((B.x.true + container.ecm$B.lx.true)/(1-container.ecm$phi))
container.ldv$centered.lrm <- container.ldv$LR.x - ((B.x.true + container.ldv$B.lx.true)/(1-container.ldv$phi))
p1 <- ggplot(subset(container.ardl, e.var == 1), aes(x = as.factor(B.lx.true), y = centered.lrm, fill = as.factor(Time) )) + geom_boxplot() + facet_wrap(~phi.label, scale = 'free') + labs(fill = "Time") + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0, color = "black", size = 1) + ylab("Estimate (centered)") + xlab("Coefficient on lag of X") + ggtitle("ARDL")
p2 <- ggplot(subset(container.ecm, e.var == 1), aes(x = as.factor(B.lx.true), y = centered.lrm, fill = as.factor(Time) )) + geom_boxplot() + facet_wrap(~phi.label, scale = 'free') + labs(fill = "Time") + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0, color = "black", size = 1) + ylab("Estimate (centered)") + xlab("Coefficient on lag of X") + ggtitle("ECM")
p3 <- ggplot(subset(container.ldv, e.var == 1), aes(x = as.factor(B.lx.true), y = centered.lrm, fill = as.factor(Time) )) + geom_boxplot() + facet_wrap(~phi.label, scale = 'free') + labs(fill = "Time") + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0, color = "black", size = 1) + ylab("Estimate (centered)") + xlab("Coefficient on lag of X") + ggtitle("LDV")
pdf("yi0-xi0-stationary-lrboxplot.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, p3, ncol = 2)
dev.off()

# combine all model datasets since geom_bar doesn't work on separate datasets
container.collapse.ardl$model.type <- "ARDL"
container.collapse.ecm$model.type <- "ECM"
container.collapse.ldv$model.type <- "LDV"
container.collapse.static$model.type <- "Static"

# Create two appended datasets: 1 short run, 1 long run:
# short run: time, e.var, phi, B.lx.true, reject.beta.x, reject.beta.x.truth, mse.beta, model.type
container.collapse.all.SR <- rbind(container.collapse.ardl[c(1,2,3,4,5,7,9,11)], container.collapse.ecm[c(1,2,3,4,5,7,9,11)], container.collapse.ldv[c(1,2,3,4,5,7,9,11)], container.collapse.static[c(1,2,3,4,5,6,7,8)])

# long run: time, e.var, phi, B.lx.true, reject.lrm, reject.lrm.truth, median square error, model type
container.collapse.all.LR <- rbind(container.collapse.ardl[c(1,2,3,4,6,8,10,11)], container.collapse.ecm[c(1,2,3,4,6,8,10,11)], container.collapse.ldv[c(1,2,3,4,6,8,10,11)])
container.collapse.all.LR$Time <- as.factor(container.collapse.all.LR$Time)
container.collapse.all.SR$Time <- as.factor(container.collapse.all.SR$Time)
container.collapse.all.LR$e.var <- as.factor(container.collapse.all.LR$e.var)
container.collapse.all.SR$e.var <- as.factor(container.collapse.all.SR$e.var)
container.collapse.all.SR$phi <- as.factor(container.collapse.all.SR$phi)
container.collapse.all.LR$phi <- as.factor(container.collapse.all.LR$phi)
container.collapse.all.SR$B.lx.true <- as.factor(container.collapse.all.SR$B.lx.true)
container.collapse.all.LR$B.lx.true <- as.factor(container.collapse.all.LR$B.lx.true)
container.collapse.all.SR$mc.type[container.collapse.all.SR$phi == .8 & container.collapse.all.SR$B.lx.true == 1] <- "Alpha = 0.8, Lag X effect = 1"
container.collapse.all.SR$mc.type[container.collapse.all.SR$phi == .8 & container.collapse.all.SR$B.lx.true == -1] <- "Alpha = 0.8, Lag X effect = -1"
container.collapse.all.SR$mc.type[container.collapse.all.SR$phi == .2 & container.collapse.all.SR$B.lx.true == 1] <- "Alpha = 0.2, Lag X effect = 1"
container.collapse.all.SR$mc.type[container.collapse.all.SR$phi == .2 & container.collapse.all.SR$B.lx.true == -1] <- "Alpha = 0.2, Lag X effect = -1"
table(container.collapse.all.SR$mc.type)
container.collapse.all.LR$mc.type[container.collapse.all.LR$phi == .8 & container.collapse.all.LR$B.lx.true == 1] <- "Alpha = 0.8, Lag X effect = 1"
container.collapse.all.LR$mc.type[container.collapse.all.LR$phi == .8 & container.collapse.all.LR$B.lx.true == -1] <- "Alpha = 0.8, Lag X effect = -1"
container.collapse.all.LR$mc.type[container.collapse.all.LR$phi == .2 & container.collapse.all.LR$B.lx.true == 1] <- "Alpha = 0.2, Lag X effect = 1"
container.collapse.all.LR$mc.type[container.collapse.all.LR$phi == .2 & container.collapse.all.LR$B.lx.true == -1] <- "Alpha = 0.2, Lag X effect = -1"
table(container.collapse.all.LR$mc.type)

# Rejection of True estimates: ------------
# now create a geom_bar when error variance = 1:
p1 <- ggplot(data = subset(container.collapse.all.SR, e.var == 1), aes(x = model.type, y = reject.beta.x.truth, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0.05, color = "black", size = 1)  + ylab("Proportion Rejected") + xlab("Model") + ggtitle("Short-Run") + facet_wrap(~mc.type) + scale_y_continuous(limits = c(0, 0.1))

p2 <- ggplot(data = subset(container.collapse.all.LR, e.var == 1), aes(x = model.type, y = reject.lrm.truth, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0.05, color = "black", size = 1) + ylab("Proportion Rejected") + xlab("Model") + ggtitle("Long-Run") + facet_wrap(~mc.type)

pdf("yi0-xi0-stationary-rejecttruth.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, ncol = 2)
dev.off()

# now create a geom_bar when error variance = 5:
p1 <- ggplot(data = subset(container.collapse.all.SR,  e.var == 5), aes(x = model.type, y = reject.beta.x.truth, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0.05, color = "black", size = 1)  + ylab("Proportion Rejected") + xlab("Model") + ggtitle("Short-Run") + facet_wrap(~mc.type) + scale_y_continuous(limits = c(0, 0.1))

p2 <- ggplot(data = subset(container.collapse.all.LR, e.var == 5), aes(x = model.type, y = reject.lrm.truth, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0.05, color = "black", size = 1) + ylab("Proportion Rejected") + xlab("Model") + ggtitle("Long-Run") + facet_wrap(~mc.type)

pdf("yi0-xi0-stationary-rejecttruth-evar5.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, ncol = 2)
dev.off()

# --- Rejection of = 0 (type II error) -------------
# evar = 1
p1 <- ggplot(data = subset(container.collapse.all.SR, e.var == 1), aes(x = model.type, y = reject.beta.x, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0.95, color = "black", size = 1)  + ylab("Proportion Rejected") + xlab("Model") + ggtitle("Short-Run") + facet_wrap(~mc.type)

p2 <- ggplot(data = subset(container.collapse.all.LR, e.var == 1), aes(x = model.type, y = reject.lrm, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0.95, color = "black", size = 1) + ylab("Proportion Rejected") + xlab("Model") + ggtitle("Long-Run") + facet_wrap(~mc.type)

pdf("yi0-xi0-stationary-reject0.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, ncol = 2)
dev.off()

# evar = 5
p1 <- ggplot(data = subset(container.collapse.all.SR,  e.var == 5), aes(x = model.type, y = reject.beta.x, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0.95, color = "black", size = 1)  + ylab("Proportion Rejected") + xlab("Model") + ggtitle("Short-Run") + facet_wrap(~mc.type) 

p2 <- ggplot(data = subset(container.collapse.all.LR, e.var == 5), aes(x = model.type, y = reject.lrm, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + geom_hline(yintercept = 0.95, color = "black", size = 1) + ylab("Proportion Rejected") + xlab("Model") + ggtitle("Long-Run") + facet_wrap(~mc.type)

pdf("yi0-xi0-stationary-reject0-evar5.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, ncol = 2)
dev.off()

# ----- Now MSE across e.var and time, for SR: -------
# evar = 1, short run
p1 <- ggplot(data = subset(container.collapse.all.SR,  e.var == 1), aes(x = model.type, y = mse.beta, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + ylab("MSE") + xlab("Model") + ggtitle("Variance = 1") + facet_wrap(~mc.type) 

# evar = 5, short-run
p2 <- ggplot(data = subset(container.collapse.all.SR,  e.var == 5), aes(x = model.type, y = mse.beta, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + ylab("MSE") + xlab("Model") + ggtitle("Variance = 5") + facet_wrap(~mc.type) 

pdf("yi0-xi0-stationary-mse-sr.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, ncol = 2)
dev.off()

# ------ Median square error across e.var and time, for LR: ---------
# evar = 1, long run
p1 <- ggplot(data = subset(container.collapse.all.LR,  e.var == 1), aes(x = model.type, y = mse.lrm, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + ylab("Median Square Error") + xlab("Model") + ggtitle("Variance = 1") + facet_wrap(~mc.type)

# evar = 5, long-run
p2 <- ggplot(data = subset(container.collapse.all.LR,  e.var == 5), aes(x = model.type, y = mse.lrm, fill = Time)) + geom_bar(stat = "identity", position = position_dodge()) + theme_minimal() + scale_fill_brewer(palette="Oranges") + ylab("Median Square Error") + xlab("Model") + ggtitle("Variance = 5") + facet_wrap(~mc.type)

pdf("yi0-xi0-stationary-mse-lr.pdf", width = 1.618*7, height = 7)
grid.arrange(p1, p2, ncol = 2)
dev.off()

# double check that EcM/ARDL rsults are identical:
cbind(container.collapse.ardl$mse.lrm,container.collapse.ecm$mse.lrm) # yup